C     ---------------------------------------------------------------
C
C     Postprocessing: calculate various cross sections
C
C     computing total and differential extinction and absorption
C     cross-sections, using the formulas given by Mackowski, Proc. R.
C     Soc. Lond. A 433, 599 (1991) and Xu, Appl. Opt. 36, 9496 (1997)
C
C     ---------------------------------------------------------------
      subroutine cross_section_dipole(nL, cext, cabs, csca, cscd, ek,
     $     drot,nmax,uvmax, k,r0,fint, atr1, btr1, as,bs, pp0, qq0)      
      implicit double precision(a-h, o-z)
      include 'gmm01f.par'
      parameter (nmp=np*(np+2),nmp0=(np+1)*(np+4)/2)
      parameter (ni0=np*(np+1)*(2*np+1)/3+np*np)
      parameter (nrc=4*np*(np+1)*(np+2)/3+np)
      parameter (nij=nLp*(nLp-1)/2)

      integer u, u0, nmax(nLp), uvmax(nLp), ind(nLp)
      double precision k,r0(6,nLp),drot(nrc,nij),csca,cext,cabs,cscd
      double precision cscai(nLp), cexti(nLp),cabsi(nLp), cscdi(nLp)
      double precision rsr(np,nLp),rsi(np,nLp),rsx(np,nLp),px(np,nLp)
      complex*16 atr1(ni0,nij),btr1(ni0,nij),at(nmp),bt(nmp),ek(np,nij)
      complex*16 as(nLp,nmp),bs(nLp,nmp),as1(nLp,nmp),bs1(nLp,nmp)
      complex*16 A,B,A0,B0, pp0(nLp,nmp),qq0(nLp,nmp)

C     constants
      pih    = dacos(0.d0)
      twopi  = 4.d0*pih
      pione  = 2.d0*pih

C     初始化需要计算的物理量
      cext = 0.d0
      cabs = 0.d0
      csca = 0.d0
C     scattering cross section with the dipole
      cscd = 0.d0
      do i = 1,nL
         cexti(i) = 0.d0
         cabsi(i) = 0.d0
         cscai(i) = 0.d0
         cscdi(i) = 0.d0
      enddo

c     the scattering cross section
      call trans(nL, r0, nmax, uvmax, fint, atr1, btr1, ek,
     +              drot, as,bs, as1,bs1, ind)

      do i=1,nL
         do imn=1,uvmax(i)
            at(imn) = as(i,imn) + as1(i,imn)
            bt(imn) = bs(i,imn) + bs1(i,imn)
         enddo
         do n=1,nmax(i)
            sc = 0.d0
            iL = n * n
            do m=-n,n
               sc =  sc + dble(dconjg(as(i,iL))*at(iL))
               sc =  sc + dble(dconjg(bs(i,iL))*bt(iL))
               iL = iL + 1
            enddo
            cscdi(i) = cscdi(i) + sc
         enddo
      enddo

c we calculate the scattered power (excluding the dipole)
      nL = nL - 1
      call trans(nL,r0,nmax,uvmax,fint,atr1,btr1,ek,
     +              drot,as,bs,as1,bs1,ind)
      do i=1,nL
         do imn=1,uvmax(i)
            at(imn)=as(i,imn)+as1(i,imn)
            bt(imn)=bs(i,imn)+bs1(i,imn)
         enddo
         do n=1,nmax(i)
            sc = 0.d0
            iL = n * n
            do m=-n,n
               sc = sc + dble(dconjg(as(i,iL))*at(iL))
               sc = sc + dble(dconjg(bs(i,iL))*bt(iL))
               iL = iL + 1
            enddo
            cscai(i) = cscai(i) + sc
         enddo
      enddo

c  ---------------------------------------------------------------
c  computing total and differential extinction and absorption
c  cross-sections, using the formulas given by Mackowski, Proc. R.
c  Soc. Lond. A 433, 599 (1991) and Xu, Appl. Opt. 36, 9496 (1997)
c  ---------------------------------------------------------------
c we calculate the power absorbed or scattered from the incident wave
c using eq. 47 of the manual
      do i=1,nL
         do imn=1,uvmax(i)
c            write (6, *) i, imn, pp0(i,imn), qq0(i,imn)
            A = DCONJG(pp0(i,imn)) * as(i,imn)
            B = DCONJG(qq0(i,imn)) * bs(i,imn)
            cexti(i) = cexti(i) + dble(A + B)
         enddo
         cabsi(i) = cexti(i) - cscai(i)
      enddo
C
      nL = nL + 1
      cscai(nL) = 0.0D0
      cexti(nL) = 0.0D0

c     sum individual particle contributions
c     and convert to relative rates normalized to the power emitted
c     by a free dipole:
c        The power emitted by a free dipole is
c        P0 = 4/3 pi |E0|^2/(k mu0 omega)
c        The power emitted by the dipole in the aggregate is given by
c        (xuy95, Gl. 57)
c        P = 4 pi/k^2 sigma = (k|E0|^2)/(2 omega mu0) * 4 pi/k^2 * sum
c        The rates are given by ri = Pi/(hbar omega)
c        therefore r/r0 = P/P0 = 3/2 sum
      cscd = 0.0D0
      csca = 0.0D0
      cext = 0.0D0
      do i=1,nL
         cscdi(i) = 3.0d0/2.0d0 * cscdi(i)
         cscai(i) = 3.0d0/2.0d0 * cscai(i)
         cexti(i) = 3.0d0/2.0d0 * cexti(i)
         cabsi(i) = 3.0d0/2.0d0 * cabsi(i)
         cscd = cscd + cscdi(i)
         csca = csca + cscai(i)
         cext = cext + cexti(i)
      enddo
      cabs = cext - csca

      return 
      end
